MODULE optimize_contractions
  !---------------------------------------------------------------------
  ! Completeness-optimized contraction pattern generation procedures
  !---------------------------------------------------------------------
  ! Copyright (C) 2004-2008 Pekka Manninen
  !                    2010 Jussi Lehtola
  !
  ! This program is distributed under the terms of the GNU General
  ! Public License as published by the Free Software Foundation;
  ! either version 2 of the License, or (at your option) any later version.
  !
  ! This program is distributed in the hope that it will be useful,
  ! but WITHOUT ANY WARRANTY; without even the implied warranty of
  ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. In addition, the 
  ! author must be properly cited whenever the program or some parts 
  ! originated on it are used. 
  !
  ! Please note, that the author must be properly cited whenever 
  ! this program or some parts originated on it are employed. 
  !
  !---------------------------------------------------------------------
  USE definitions
  USE error
  USE random
  USE completeness_profile, ONLY : comp_prof
  USE measure
  USE textutil

  IMPLICIT NONE

  PRIVATE
  PUBLIC :: opt_contr, opt_contr_simann, opt_contr_lottery

CONTAINS

  !-------------------------------------------------------------------  
  ! Routines for contraction pattern optimization
  !-------------------------------------------------------------------

  SUBROUTINE opt_contr(l, optmin, optmax, X, C, usesto, measure, conv)
    ! Optimize the contractions

    IMPLICIT NONE
    ! Value of angular momentum
    INTEGER, INTENT(IN) :: l
    ! Exponent range to optimize in
    REAL(KIND=prec), INTENT(IN) :: optmin, optmax
    ! Exponents (hopefully completeness-optimized)
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: X
    ! Contraction coefficients to optimize
    REAL(KIND=prec), DIMENSION(:,:), INTENT(INOUT) :: C
    ! Use STOs to probe?
    LOGICAL, INTENT(IN) :: usesto
    ! Which measure to use
    INTEGER, INTENT(IN) :: measure
    ! The given convergence threshold
    REAL(KIND=prec), INTENT(IN) :: conv

    ! Probing exponents
    REAL(KIND=prec), DIMENSION(optpoints) :: alphas

    ! Evolution populations
    REAL(KIND=prec), DIMENSION(SIZE(C,1),SIZE(C,2),pop_dim*SIZE(C,1)*SIZE(C,2)) :: popul
    ! Values of cost function for the populations
    REAL(KIND=prec), DIMENSION(pop_dim*SIZE(C,1)*SIZE(C,2)) :: cost

    ! Time variables
    INTEGER :: time1, time2, cpucount

    ! Amount of exponents and contractions
    INTEGER :: nx, nc

    ! Location of minimum cost
    INTEGER, DIMENSION(1) :: minind

    nx=SIZE(X)
    nc=SIZE(C,2)

    ! Fill probing exponents
    CALL get_alphas(optmin,optmax,alphas)

    WRITE(*,'(A,I3,A)') 'Optimizing a set of ', nc, ' contractions. This may take *a lot* of time.'

    ! Get wall clock time
    CALL system_clock(time1,cpucount)
    
    ! Make a population of contractions
    CALL twod_init(l, measure, X, popul, cost, alphas, usesto)
    ! Run the optimization
    CALL twod_opt(l, measure, X, popul, cost, alphas, usesto, conv)

    ! Get wall clock time
    CALL system_clock(time2,cpucount)

    WRITE (*,'(A,F5.1,A,/)') 'Optimization took ',(time2-time1)*1.0/cpucount,' seconds.'
    
    ! Get minimum cost value
    minind=MINLOC(cost)

    ! Get optimal set
    C(1:nx,1:nc)=popul(1:nx,1:nc,minind(1))

    ! Print values
    WRITE(*,'(A)') 'Values for optimized set:'
    CALL f_everything(l, X, C, alphas, usesto)
    
  END SUBROUTINE opt_contr


  SUBROUTINE opt_contr_simann(l, optmin, optmax, X, C, usesto, measure, Nsteps)
    ! Optimize the contractions by simulated annealing

    IMPLICIT NONE
    ! Value of angular momentum
    INTEGER, INTENT(IN) :: l
    ! Exponent range to optimize in
    REAL(KIND=prec), INTENT(IN) :: optmin, optmax
    ! Exponents (hopefully completeness-optimized)
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: X
    ! Contraction coefficients to optimize
    REAL(KIND=prec), DIMENSION(:,:), INTENT(INOUT) :: C
    ! Use STOs to probe?
    LOGICAL, INTENT(IN) :: usesto
    ! Which measure to use
    INTEGER, INTENT(IN) :: measure
    ! Number of steps in each annealing cycle
    INTEGER, INTENT(IN) :: Nsteps

    ! Random number generator, defined in gslrand.c
    REAL(KIND=prec) :: rnd_gaussian

    ! Probing exponents
    REAL(KIND=prec), DIMENSION(optpoints) :: alphas

    ! Time variables
    INTEGER :: time1, time2, cpucount

    ! Amount of exponents and contractions
    INTEGER :: nx, nc

    ! Loop indices
    INTEGER :: i, j

    nx=SIZE(X)
    nc=SIZE(C,2)

    ! Fill probing exponents
    CALL get_alphas(optmin,optmax,alphas)

    WRITE(*,'(A,I3,A)') 'Optimizing a set of ', nc, ' contractions. This may take *a lot* of time.'

    ! Get wall clock time
    CALL system_clock(time1,cpucount)

    ! Initial values for contraction
    DO j = 1, SIZE(C,2)
       DO i = 1, SIZE(C,1)
          C(i,j) = rnd_gaussian()
       END DO
    END DO

    ! Run simulated annealing
    CALL twod_simann(l, measure, X, C, alphas, usesto, Nsteps)

    ! Get wall clock time
    CALL system_clock(time2,cpucount)

    WRITE (*,'(A,F5.1,A,/)') 'Optimization took ',(time2-time1)*1.0/cpucount,' seconds.'
    
    ! Print values
    WRITE(*,'(A)') 'Values for optimized set:'
    CALL f_everything(l, X, C, alphas, usesto)
    
  END SUBROUTINE opt_contr_simann


  SUBROUTINE opt_contr_lottery(l, optmin, optmax, X, C, usesto, measure, Nsteps)
    ! Optimize the contractions by lottery

    IMPLICIT NONE
    ! Value of angular momentum
    INTEGER, INTENT(IN) :: l
    ! Exponent range to optimize in
    REAL(KIND=prec), INTENT(IN) :: optmin, optmax
    ! Exponents (hopefully completeness-optimized)
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: X
    ! Contraction coefficients to optimize
    REAL(KIND=prec), DIMENSION(:,:), INTENT(INOUT) :: C
    ! Use STOs to probe?
    LOGICAL, INTENT(IN) :: usesto
    ! Which measure to use
    INTEGER, INTENT(IN) :: measure
    ! Number of lotery phases
    INTEGER, INTENT(IN) :: Nsteps

    ! Random number generator, defined in gslrand.c
    REAL(KIND=prec) :: rnd_gaussian

    ! Probing exponents
    REAL(KIND=prec), DIMENSION(optpoints) :: alphas

    ! Time variables
    INTEGER :: time1, time2, cpucount

    ! Amount of exponents and contractions
    INTEGER :: nx, nc

    ! Loop indices
    INTEGER :: i, j

    nx=SIZE(X)
    nc=SIZE(C,2)

    ! Fill probing exponents
    CALL get_alphas(optmin,optmax,alphas)

    WRITE(*,'(A,I3,A)') 'Optimizing a set of ', nc, ' contractions. This may take *a lot* of time.'

    ! Get wall clock time
    CALL system_clock(time1,cpucount)

    ! Initial values for contraction
    DO j = 1, SIZE(C,2)
       DO i = 1, SIZE(C,1)
          C(i,j) = rnd_gaussian()
       END DO
    END DO

    ! Run simulated annealing
    CALL twod_lottery(l, measure, X, C, alphas, usesto, Nsteps, 100)

    ! Get wall clock time
    CALL system_clock(time2,cpucount)

    WRITE (*,'(A,F5.1,A,/)') 'Optimization took ',(time2-time1)*1.0/cpucount,' seconds.'

    ! Print values
    WRITE(*,'(A)') 'Values for optimized set:'
    CALL f_everything(l, X, C, alphas, usesto)
    
  END SUBROUTINE opt_contr_lottery



  SUBROUTINE twod_init(l,measure,X,pop,cost,alphas,usesto)
    IMPLICIT NONE

    ! Angular momentum type
    INTEGER, INTENT(IN) :: l
    ! Measure to use in optimization
    INTEGER, INTENT(IN) :: measure
    ! Optimal exponents
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: X
    ! Population of contraction coefficients
    REAL(KIND=prec), DIMENSION(:,:,:), INTENT(OUT) :: pop
    ! Values of cost function for population
    REAL(KIND=prec), DIMENSION(:), INTENT(OUT) :: cost
    ! Scanning exponents
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: alphas
    ! Use STO to probe?
    LOGICAL, INTENT(IN) :: usesto

    ! Gaussian random number generator, defined in gslrand.c
    REAL(KIND=prec) :: rnd_gaussian
    ! Loop indices
    INTEGER :: i, j, k

    IF(SIZE(cost) /= SIZE(pop,3)) STOP 'Cost vector of wrong size'

    ! Set contractions
    DO k = 1, SIZE(pop,3)
       DO j = 1, SIZE(pop,2)
          DO i = 1, SIZE(pop,1)
             pop(i,j,k) = rnd_gaussian()
          END DO
       END DO
    END DO

    ! Calculate costs
    SELECT CASE(measure)
    CASE(1) ! Use maximum deviance
       DO k = 1, SIZE(pop,3)
          cost(k) = f_maxdev(l,X,pop(:,:,k),alphas,usesto)
       END DO
    CASE(2) ! Use RMS deviance
       DO k = 1, SIZE(pop,3)
          cost(k) = f_lsq(l,X,pop(:,:,k),alphas,usesto)
       END DO
    CASE(3) ! Use integral
       DO k = 1, SIZE(pop,3)
          cost(k) = f_sum(l,X,pop(:,:,k),alphas,usesto)
       END DO
    CASE DEFAULT
       CALL error_stop('Error in optimize_contractions.f90 - no such measure!')
    END SELECT
  END SUBROUTINE twod_init


  SUBROUTINE twod_opt(l, measure, X, pop, cost, alphas, usesto, conv)
    ! differential evolution optimization routine for matrices
    IMPLICIT NONE

    ! Angular momentum
    INTEGER, INTENT(IN) :: l
    ! Measure to minimize
    INTEGER, INTENT(IN) :: measure
    ! Exponents
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: X
    ! The population of contraction coefficients
    REAL(KIND=prec), DIMENSION(:,:,:), INTENT(INOUT) :: pop
    ! Cost functions for the populations
    REAL(KIND=prec), DIMENSION(:), INTENT(INOUT) :: cost
    ! Scanning exponents
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: alphas
    ! Use STO to probe?
    LOGICAL, INTENT(IN) :: usesto
    ! Convergence criterion
    REAL(KIND=prec), INTENT(IN) :: conv

    ! New population
    REAL(KIND=prec), DIMENSION(SIZE(pop,1), SIZE(pop,2), SIZE(pop,3)) :: new_pop
    ! Trial coefficients
    REAL(KIND=prec), DIMENSION(SIZE(pop,1), SIZE(pop,2)) :: trial


    ! Number of exponents and contractions
    INTEGER :: nx, nc
    ! Population size
    INTEGER :: pop_size

    ! Helpers
    INTEGER :: iter, i, idx1, idx2, idx3

    ! Values of cost function
    REAL(KIND=prec) :: score, oldit

    ! RND variable for determining which exponents to pair
    REAL(KIND=prec), DIMENSION(SIZE(pop,1),SIZE(pop,2)) :: p

    ! Has the evolution converged?
    LOGICAL :: converged

    nx = SIZE(pop,1)
    nc = SIZE(pop,2)
    pop_size = SIZE(pop,3)

    IF(SIZE(cost) /= pop_size) THEN
       STOP 'Vector cost is of wrong size (twod_opt).'
    END IF
    
    new_pop = 0.0_prec
    oldit = 1.0_prec
    score=0.0_prec

    SELECT CASE(measure)
    CASE(1)
       WRITE(*,'(/,A10,A10,A15)') ' ', 'Iteration', 'Max deviation'
    CASE(2)
       WRITE(*,'(/,A10,A10,A15)') ' ', 'Iteration', 'RMS deviation'
    CASE(3)
       WRITE(*,'(/,A10,A10,A15)') ' ', 'Iteration', 'Integral'
    END SELECT

    iter = 1
    converged = .FALSE.
    DO WHILE (iter <= maxit .AND. .NOT. converged)
       DO i = 1, pop_size
          CALL rnd_triple(i,pop_size,idx1,idx2,idx3)
          p = rnd_mat(nx,nc)
          WHERE (p<evol_r)
             trial = pop(:,:,idx3) + evol_c*(pop(:,:,idx1)-pop(:,:,idx2))
          ELSEWHERE
             trial = pop(:,:,i)
          END WHERE

          ! Calculate cost
          SELECT CASE(measure)
          CASE(1)
             ! Use maximum deviance
             score = f_maxdev(l,X,trial,alphas,usesto)
          CASE(2)
             ! Use squared deviance
             score = f_lsq(l,X,trial,alphas,usesto)
          CASE(3)
             ! Use integral
             score = f_sum(l,X,trial,alphas,usesto)
          CASE DEFAULT
             CALL error_stop('Error in twod_opt - invalid measure!')
          END SELECT

          IF(score < cost(i)) THEN
             new_pop(:,:,i) = trial
             cost(i) = score
          ELSE
             new_pop(:,:,i) = pop(:,:,i)
          END IF
       END DO

       ! Change populations
       pop = new_pop
       IF(MOD(iter,istep) == 0) THEN
          WRITE(*,'(A10,I10,E15.6)') ' ', iter, MINVAL(cost)

          IF (oldit - MINVAL(cost) < conv) THEN
             converged = .TRUE.
          END IF

          oldit = MINVAL(cost)
       END IF

       iter = iter + 1
    END DO

    IF (iter < maxit) THEN
       WRITE(*,'(A)') 'Optimization converged!'
    ELSE
       WRITE(*,'(A)') 'NOT CONVERGED: Maximum number of iterations reached.'
    END IF
  END SUBROUTINE twod_opt

  SUBROUTINE twod_simann(l, measure, X, init_C, alphas, usesto, Nsteps)
    ! Simulated annealing routine for optimization
    IMPLICIT NONE

    ! Angular momentum
    INTEGER, INTENT(IN) :: l
    ! Measure to minimize
    INTEGER, INTENT(IN) :: measure
    ! Exponents
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: X
    ! Contraction coefficients
    REAL(KIND=prec), DIMENSION(:,:), INTENT(INOUT) :: init_C
    ! Scanning exponents
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: alphas
    ! Use STO to probe?
    LOGICAL, INTENT(IN) :: usesto
    ! Amount of steps to run
    INTEGER, INTENT(IN) :: Nsteps

    ! Random number generator, defined in gslrand.c
    REAL(KIND=prec) :: rnd_uniform

    ! Number of exponents and contractions
    INTEGER :: nx, nc

    ! Trial coefficients
    REAL(KIND=prec), DIMENSION(SIZE(init_C,1),SIZE(init_C,2)) :: trial_C
    REAL(KIND=prec) :: trial_cost
    ! Current coefficients
    REAL(KIND=prec), DIMENSION(SIZE(init_C,1),SIZE(init_C,2)) :: cur_C
    REAL(KIND=prec) :: cur_cost
    ! Minimum value of cost and minimal coefficients
    REAL(KIND=prec), DIMENSION(SIZE(init_C,1),SIZE(init_C,2)) :: min_C
    REAL(KIND=prec) :: min_cost

    ! Cost of past iteration
    REAL(KIND=prec) :: old_cost

    ! Maximum step size
    REAL(KIND=prec), PARAMETER :: stepsize = 1
    ! Step size for trials, dependent on temperature
    REAL(KIND=prec) :: step_T    

    ! Loop variables
    INTEGER :: istep, icycle

    ! Number of temperature cycles to run
    INTEGER :: Ncycles

    ! Has the evolution converged?
    LOGICAL :: converged

    !! Temperature variables
    ! The cost function should be <= 1 for all basis sets
    ! Set initial temperature to 0.1
    REAL(KIND=prec), PARAMETER :: T_max = 0.0001_prec
    REAL(KIND=prec) :: T, T0

    nx = SIZE(init_C,1)
    nc = SIZE(init_C,2)

    T0=T_max

    ! Number of cycles is determined by how many cycles are
    ! needed to drop T0 below a reasonable limit, here 10^-10
    Ncycles=CEILING((10.0_prec+LOG10(T_max))/(-LOG10(0.8_prec))+1.0_prec)

    ! Normalize initial contraction
!    CALL normalize_contraction(init_C)

    ! Calculate initial cost value
    cur_C=init_C
    cur_cost=0.0_prec
    trial_cost=0.0_prec
    old_cost=0.0_prec
    SELECT CASE(measure)
    CASE(1)
       ! Use maximum deviance
       cur_cost = f_maxdev(l,X,cur_C,alphas,usesto)
    CASE(2)
       ! Use squared deviance
       cur_cost = f_lsq(l,X,cur_C,alphas,usesto)
    CASE(3)
       ! Use integral
       cur_cost = f_sum(l,X,cur_C,alphas,usesto)
    CASE DEFAULT
       CALL error_stop('Error in twod_simann - invalid measure!')
    END SELECT

    ! Set initial values

    min_cost=cur_cost
    min_C=cur_C

    SELECT CASE(measure)
    CASE(1)
       WRITE(*,'(/,A10,A10,A15)') ' ', 'Iteration', 'Max deviation'
    CASE(2)
       WRITE(*,'(/,A10,A10,A15)') ' ', 'Iteration', 'RMS deviation'
    CASE(3)
       WRITE(*,'(/,A10,A10,A15)') ' ', 'Iteration', 'Integral'
    END SELECT

    converged = .FALSE.
    ! Loop over annealing cycles
    DO icycle=1,Ncycles
       ! Loop over single cycle
       DO istep=1,Nsteps

          ! Calculate temperature
          T=0.1_prec*T0 + 0.9_prec*T0*(istep-1.0_prec)/(Nsteps-1.0_prec)

          ! Calculate step size
          step_T=T*stepsize

          ! Make trial
          trial_C=cur_C+step_T*(2.0_prec*rnd_mat(nx,nc)-1)

          ! Normalize trial
!          CALL normalize_contraction(trial_C)

          ! Calculate cost of trial
          SELECT CASE(measure)
          CASE(1)
             ! Use maximum deviance
             trial_cost = f_maxdev(l,X,trial_C,alphas,usesto)
          CASE(2)
             ! Use squared deviance
             trial_cost = f_lsq(l,X,trial_C,alphas,usesto)
          CASE(3)
             ! Use integral
             trial_cost = f_sum(l,X,trial_C,alphas,usesto)
          CASE DEFAULT
             CALL error_stop('Error in twod_simann - invalid measure!')
          END SELECT

          ! Determine whether to accept trial
          IF(trial_cost<cur_cost) THEN
             ! Smaller cost, accept automatically
             cur_C=trial_C
             cur_cost=trial_cost
             ! Did we go under the minimum value reached so far?
             IF(trial_cost<min_cost) THEN
                min_C=trial_C
                min_cost=trial_cost
             END IF
          ELSEIF (rnd_uniform() < exp((trial_cost-cur_cost)/T)) THEN
             ! Higher cost, but accepted Metropolis trial
             cur_C=trial_C
             cur_cost=trial_cost
          END IF
       END DO

       ! Print out current status
       IF(measure==2) THEN
          WRITE(*,'(A10,I10,E15.6,E15.6,E15.6,E15.6)') ' ', icycle, &
               SQRT(cur_cost), SQRT(cur_cost)-SQRT(old_cost), T0, SQRT(min_cost)
       ELSE
          WRITE(*,'(A10,I10,E15.6,E15.6,E15.6,E15.6)') ' ', icycle, &
               ABS(cur_cost), ABS(cur_cost)-ABS(old_cost), T0, ABS(min_cost)
       END IF
       
       ! Lower reference
       T0=0.8_prec*T0
       
       ! Store value of old iteration
       old_cost=cur_cost
    END DO

    ! Set return coefficients to minimal ones
    init_C=min_C

  END SUBROUTINE twod_simann

  SUBROUTINE twod_lottery(l, measure, X, init_C, alphas, usesto, Nsteps, Ncycles)
    ! Lottery optimization
    IMPLICIT NONE

    ! Angular momentum
    INTEGER, INTENT(IN) :: l
    ! Measure to minimize
    INTEGER, INTENT(IN) :: measure
    ! Exponents
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: X
    ! Contraction coefficients
    REAL(KIND=prec), DIMENSION(:,:), INTENT(INOUT) :: init_C
    ! Scanning exponents
    REAL(KIND=prec), DIMENSION(:), INTENT(IN) :: alphas
    ! Use STO to probe?
    LOGICAL, INTENT(IN) :: usesto
    ! Amount of steps to run
    INTEGER, INTENT(IN) :: Nsteps
    ! Number of temperature cycles to run
    INTEGER, INTENT(IN) :: Ncycles

    ! Random number generator, defined in gslrand.c
    REAL(KIND=prec) :: rnd_gaussian

    ! Number of exponents and contractions
    INTEGER :: nx, nc

    ! Trial coefficients
    REAL(KIND=prec), DIMENSION(SIZE(init_C,1),SIZE(init_C,2)) :: trial_C
    REAL(KIND=prec) :: trial_cost=0.0_prec
    ! Minimum value of cost
    REAL(KIND=prec) :: min_cost=0.0_prec

    ! Cost of past iteration
    REAL(KIND=prec) :: old_cost=1.0_prec

    ! Loop variables
    INTEGER :: istep, icycle, i, j


    nx = SIZE(init_C,1)
    nc = SIZE(init_C,2)

    ! Get initial cost
    SELECT CASE(measure)
    CASE(1)
       ! Use maximum deviance
       min_cost = f_maxdev(l,X,init_C,alphas,usesto)
    CASE(2)
       ! Use squared deviance
       min_cost = f_lsq(l,X,init_C,alphas,usesto)
    CASE(3)
       ! Use integral
       min_cost = f_sum(l,X,init_C,alphas,usesto)
    CASE DEFAULT
       CALL error_stop('Error in twod_lottery - invalid measure!')
    END SELECT

    ! Print out
    SELECT CASE(measure)
    CASE(1)
       WRITE(*,'(/,A10,A10,A15)') ' ', 'Iteration', 'Max deviation'
    CASE(2)
       WRITE(*,'(/,A10,A10,A15)') ' ', 'Iteration', 'RMS deviation'
    CASE(3)
       WRITE(*,'(/,A10,A10,A15)') ' ', 'Iteration', 'Integral'
    END SELECT

    
    DO icycle=1,Ncycles
       DO istep=1,Nsteps
          
          ! Generate a trial
          DO j = 1, nc
             DO i = 1, nx
                trial_C(i,j) = rnd_gaussian()
             END DO
          END DO

          ! and calculate its cost
          SELECT CASE(measure)
          CASE(1)
             ! Use maximum deviance
             trial_cost = f_maxdev(l,X,trial_C,alphas,usesto)
          CASE(2)
             ! Use squared deviance
             trial_cost = f_lsq(l,X,trial_C,alphas,usesto)
          CASE(3)
             ! Use integral
             trial_cost = f_sum(l,X,trial_C,alphas,usesto)
          CASE DEFAULT
             CALL error_stop('Error in twod_lottery - invalid measure!')
          END SELECT

          ! If cost is smaller than minimal cost, update values
          IF(trial_cost < min_cost) THEN
             init_C = trial_C
             min_cost = trial_cost
          END IF
       END DO

       ! Print out current status
       IF(measure==2) THEN
          WRITE(*,'(A10,I10,E15.6,E15.6)') ' ', icycle, &
               SQRT(min_cost), SQRT(min_cost)-SQRT(min_cost)
       ELSE
          WRITE(*,'(A10,I10,E15.6,E15.6)') ' ', icycle, &
               ABS(min_cost), ABS(min_cost)-ABS(min_cost)
       END IF

       ! Store value of old iteration
       old_cost=min_cost
    END DO

  END SUBROUTINE twod_lottery

END MODULE optimize_contractions
